
#include <stdio.h>
#include <math.h>
#include "../lib/inc/signal_process.h"

/**
*/

int main()
{
	int i, j, n;
	double a1, a2, c, c1, c2, d1, d2, q1, q2, w, w1, w2;
	double x[32], y[32], a[32], b[32];
	n = 32;
	a1 = 0.9;
	a2 = 0.3;
	x[0] = 1.0;
	y[0] = 0.0;
	for (i = 1; i < n; i++)
	{
		x[i] = a1 * x[i - 1] - a2 * y[i - 1];
		y[i] = a2 * x[i - 1] + a1 * y[i - 1];
	}
	printf("\n The Original Sequence \n");
	for (i = 0; i < n / 2; i++)
	{
		for (j = 0; j < 2; j++)
		{
			printf("	%10.7lf + j %10.7lf", x[2 * i + j], y[2 * i + j]);
		}
		printf("\n");
	}

	q1 = x[n - 1];
	q2 = y[n - 1];
	dft(x,y,a,b,n,1);
	printf("\n Discrete Fourier Transform \n");
	for (i = 0; i < n / 2; i++)
	{
		for (j = 0; j < 2; j++)
		{
			printf("	%10.7lf + j %10.7lf", a[2 * i + j], b[2 * i + j]);
		}
		printf("\n");
	}

	for (i = 0; i < n; i++)
	{
		w = 6.28318530718 / n * i;
		w1 = cos(w);
		w2 = -sin(w);
		c1 = 1.0 - a1 * w1 + a2 * w2;
		c2 = a1 * w2 + a2 * w1;
		c = c1 * c1 + c2 * c2;
		d1 = 1.0 - a1 * q1 + a2 * q2;
		d2 = a1 * q2 + a2 * q1;
		x[i] = (c1 * d1 + c2 * d2) / c;
		y[i] = (c2 * d1 - c1 * d2) / c;
	}

	printf("\n Theoretical Discrete Fourier Transform \n");
	for (i = 0; i < n / 2; i++)
	{
		for (j = 0; j < 2; j++)
		{
			printf("	%10.7lf + j %10.7lf", x[2 * i + j], y[2 * i + j]);
		}
		printf("\n");
	}
	dft(a, b, x, y, n, -1);
	printf("\n Inverse Discrete Fourier Transform \n");
	for (i = 0; i < n / 2; i++)
	{
		for (j = 0; j < 2; j++)
		{
			printf("	%10.7lf + j %10.7lf", x[2 * i + j], y[2 * i + j]);
		}
		printf("\n");
	}
	getchar();
}